home *** CD-ROM | disk | FTP | other *** search
- c imsl routine name - vbla=ddot vbdd0010
- c
- c-----------------------------------------------------------------------
- c
- c computer - vax/double
- c
- c latest revision - january 1, 1978
- c
- c purpose - compute double precision dot product
- c
- c usage - function ddot (n,dx,incx,dy,incy)
- c
- c arguments ddot - double precision sum from i=1 to n of
- c x(i)*y(i). (output)
- c x(i) and y(i) refer to specific elements
- c of dx and dy, respectively. see incx and
- c incy argument descriptions.
- c n - length of vectors x and y. (input)
- c dx - double precision vector of length
- c max(n*iabs(incx),1). (input)
- c incx - displacement between elements of dx. (input)
- c x(i) is defined to be..
- c dx(1+(i-1)*incx) if incx.ge.0 or
- c dx(1+(i-n)*incx) if incx.lt.0.
- c dy - double precision vector of length
- c max(n*iabs(incy),1). (input)
- c incy - displacement between elements of dy. (input)
- c y(i) is defined to be..
- c dy(1+(i-1)*incy) if incy.ge.0 or
- c dy(1+(i-n)*incy) if incy.lt.0.
- c
- c precision/hardware - double/all
- c
- c reqd. imsl routines - none required
- c
- c notation - information on special notation and
- c conventions is available in the manual
- c introduction or through imsl routine uhelp
- c
- c copyright - 1978 by imsl, inc. all rights reserved.
- c
- c warranty - imsl warrants only that imsl testing has been
- c applied to this code. no other warranty,
- c expressed or implied, is applicable.
- c
- c-----------------------------------------------------------------------
- c
- double precision function ddot (n,dx,incx,dy,incy)
- c
- c specifications for arguments
- double precision dx(1),dy(1)
- integer n,incx,incy
- c specifications for local variables
- integer i,m,mp1,ns,ix,iy
- c first executable statement
- ddot = 0.d0
- if (n.le.0) return
- if (incx.eq.incy) if (incx-1) 5,15,35
- 5 continue
- c code for unequal or nonpositive
- c increments.
- ix = 1
- iy = 1
- if (incx.lt.0) ix = (-n+1)*incx+1
- if (incy.lt.0) iy = (-n+1)*incy+1
- do 10 i=1,n
- ddot = ddot+dx(ix)*dy(iy)
- ix = ix+incx
- iy = iy+incy
- 10 continue
- return
- c code for both increments equal to 1.
- c clean-up loop so remaining vector
- c length is a multiple of 5.
- 15 m = n-(n/5)*5
- if (m.eq.0) go to 25
- do 20 i=1,m
- ddot = ddot+dx(i)*dy(i)
- 20 continue
- if (n.lt.5) return
- 25 mp1 = m+1
- do 30 i=mp1,n,5
- ddot = ddot+dx(i)*dy(i)+dx(i+1)*dy(i+1)+dx(i+2)*dy(i+2)+dx(i
- 1 +3)*dy(i+3)+dx(i+4)*dy(i+4)
- 30 continue
- return
- c code for positive equal increments
- c .ne.1.
- 35 continue
- ns = n*incx
- do 40 i=1,ns,incx
- ddot = ddot+dx(i)*dy(i)
- 40 continue
- return
- end
-